#include  "/home/yukp/fit_ll_angle/Event.h"
#include "/home/yukp/headfiles/Omegac2OmegamPi/RooAngularPdf.h"
vector <double> vtheta0_mc,  vtheta1_mc,  vtheta2_mc;
vector <double>   vphi0_mc,    vphi1_mc,    vphi2_mc;
vector <double> vtheta0_sig, vtheta1_sig, vtheta2_sig;
vector <double>   vphi0_sig,   vphi1_sig,   vphi2_sig;
vector <double> vweight_mc;

void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
double func0(
        vector<double> costheta, vector<double> sintheta,
        vector<double> cosphi,   vector<double> sinphi,
        double *par);
double getsigma_phsp(double *par);
Event readData(double alpha,double Phi,double alpha1,double Phi1,double alpha2,double Px,double Py,double Pz,int total_events,long SEED);
Event evts_phsp, evts_sig;
double cal_beta(double alpha,double phi);
double cal_gamma(double alpha,double phi);
TTree *save_result(TTree *t_data,double *par,double *err,double *betaerr,double *gammaerr,int flag);

/*serial number of each parameters*/
int alpha_num=0;
int phi_num=1;
int alpha1_num=2;
int phi1_num=3;
int alpha2_num=4;
int Px_num=5;
int Py_num=6;
int Pz_num=7;

double alpha_Omegac=-0.64;
double alpha_Omegam=-0.0157;
double alpha_Lambda=-0.732;
double PhiOmegac_initial=-2.1;
double PhiOmegam_initial=-2.1;

double Px_initial=0.0;
double Py_initial=0.01;
double Pz_initial=0.1;
double gen_settings[]={alpha_Omegac,PhiOmegac_initial,alpha_Omegam,PhiOmegam_initial,alpha_Lambda,Px_initial,Py_initial,Pz_initial};

/*****************************************************************************
 * end define*
 *****************************************************************************/

using namespace RooFit;
void pseudo(){

long Seed=SEEDRAND;
int total=4*30000;
{//generating the refering dataset
    cout<< "--------- generating reference toy --------- " << endl;

    evts_phsp = readData(alpha_Omegac,PhiOmegac_initial,alpha_Omegam,PhiOmegam_initial,alpha_Lambda,Px_initial,Py_initial,Pz_initial,total,Seed);
    cout<< "--- Total Events in MC = " << evts_phsp.numEvt_total   << endl;
}


    TString filename=TString("./fit_result/fit_result")+TString("_with_SEED_")+Seed+TString(".root");
    TFile *f = new TFile(filename,"RECREATE");
    TTree *t_data = new TTree("t_data","parameters");
    for (int i=0;i<1000;i++)
    {
        Seed+=1;
            {//generating the toy dataset
            cout<< "----- generating  toy ------- " << endl;
            evts_sig  = readData(alpha_Omegac,PhiOmegac_initial,alpha_Omegam,PhiOmegam_initial,alpha_Lambda,Px_initial,Py_initial,Pz_initial,total,Seed);
            cout<< "--- Total Events in SIG = " << evts_sig.numEvt_total  << endl;
            }

            cout<< "--- start fitting--- " << endl;

            //*****************************************************************************
            //* minimize *
            //*****************************************************************************

            TMinuit *gMinuit = new TMinuit(15);  //initialize TMinuit with a maximum of 10 params
            gMinuit->SetFCN(fcn);

            Double_t arglist[10];
            Int_t ierflg = 0;

            arglist[0] = 1;
            gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);



            // Set starting values and step sizes for parameters
            gMinuit->mnparm(alpha_num,  "alpha",   alpha_Omegac ,      0.01, -1, 1, ierflg);
            gMinuit->mnparm(phi_num,    "phi",     PhiOmegac_initial , 0.01, -180., 180., ierflg);
            gMinuit->mnparm(alpha1_num, "alpha1",  alpha_Omegam  ,     0.01, -1, 1, ierflg);
            gMinuit->mnparm(phi1_num,   "phi1",    PhiOmegam_initial , 0.01, -180., 180, ierflg);
            gMinuit->mnparm(alpha2_num, "alpha2",  alpha_Lambda ,   0.01, -1, 1, ierflg);
            gMinuit->mnparm(Px_num,     "Px",      Px_initial ,     0.01, -10, 10, ierflg);
            gMinuit->mnparm(Py_num,     "Py",      Py_initial ,     0.01, -10, 10, ierflg);
            gMinuit->mnparm(Pz_num,     "Pz",      Pz_initial ,     0.01, -10, 10, ierflg);

            gMinuit->FixParameter(alpha1_num);
            gMinuit->FixParameter(alpha2_num);

            // Now ready for minimization step
            arglist[0] = 500;
            arglist[1] = 1.;
            gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
            gMinuit->mnexcm("HESSE",  arglist ,1,ierflg);
            //gMinuit->mnexcm("MINOS", arglist ,1,ierflg);

            // Print results
            Double_t amin,edm,errdef;
            Int_t nvpar,nparx,icstat;
            gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);

        double par[10], err[10];
            for(int i=0; i<8; i++){
                gMinuit->GetParameter(i, par[i], err[i]);
            }


            Double_t cov_matrix[6][6];
            gMinuit->mnemat(&cov_matrix[0][0],6);

            /*0         alpha
            *1         phi
            *(2) (alpha1)
            *(3)2      phi1
            *(4) (alpha2)
            *(5)3      Px
            *(6)4      Py
            *(7)5      Pz
            */


            //calculate the progated error of beta & gamma
            double betaerr[2],gammaerr[2];
            int  position[2][2]= {{alpha_num,phi_num}, {alpha1_num,phi1_num}};
            double covariant[]= {cov_matrix[0][1],0};
            for (int i=0;i<2;i++){
                double alpha = par[position[i][0] ];
                double phi =par[ position[i][1] ];
                double alphaerr = err[position[i][0] ];
                double phierr =err[ position[i][1] ];
                double beta_derive_alpha=0;
                double gamma_derive_alpha=0;
                double beta_derive_phi=0;
                double gamma_derive_phi=0;
                if (alphaerr!=0) {
                    beta_derive_alpha=(cal_beta(alpha+alphaerr,phi)- cal_beta(alpha-alphaerr,phi))/2/alphaerr;
                    gamma_derive_alpha=(cal_gamma(alpha+alphaerr,phi)- cal_gamma(alpha-alphaerr,phi))/2/alphaerr;
                }
                if (phierr!=0) {
                    beta_derive_phi=(cal_beta(alpha,phi+phierr)- cal_beta(alpha,phi-phierr))/2/phierr;
                    gamma_derive_phi=(cal_gamma(alpha,phi+phierr)- cal_gamma(alpha,phi-phierr))/2/phierr;
                }
                betaerr[i]=sqrt(
                        pow(beta_derive_alpha,2)*pow(alphaerr,2)+
                        pow(beta_derive_phi,2)*pow(phierr,2)+
                        2*covariant[i]*beta_derive_alpha*beta_derive_phi
                           );

                gammaerr[i]=sqrt(
                        pow(gamma_derive_alpha,2)*pow(alphaerr,2)
                        +pow(gamma_derive_phi,2)*pow(phierr,2)
                        +2*covariant[i]*gamma_derive_alpha*gamma_derive_phi
                        );

                    /*    cout << "alphaerr" <<i<<" = "<<alphaerr<<endl;
                         cout << "phierr" <<i<<" = "<<phierr<<endl;
                        cout <<" beta_"<< i<<" = "<<betaerr[i]<<endl;
                        cout <<" gamma_"<< i<<" = "<<gammaerr[i]<<endl;
                */
		}
	    cout<< endl ;
cout<<i<<"st  experiment completed "<<endl;
cout<< endl ;

    t_data=save_result(t_data,par,err,betaerr,gammaerr,i);

    }
    f->Write();
    f->Close();


}





void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {


    double likelihood_0[3] = {0.,0.,0.};

    double int_sigma = getsigma_phsp(par);

    for (int i = 0; i< evts_sig.numEvt_total; i++){
        double li = -1.*log( func0(evts_sig.costhetas[i], evts_sig.sinthetas[i], evts_sig.cosphis[i], evts_sig.sinphis[i], par) );
        likelihood_0[0] += li;
    }
    likelihood_0[0]+=log(int_sigma);

    likelihood_0[2] = likelihood_0[0];


    //*************result**************//
   // cout<<" | "<<setw(10)<<"likelihood_0[2] = "<<setw(10)<<likelihood_0[2];
    f = 2. * (likelihood_0[2]);

    /*cout<<" | "<<setw(10)<<f;
    printf("%20.10f",f);
    cout<<endl;
    */
}


double getsigma_phsp(double *par){
    double sigma = 0;
    for (int i = 0; i< evts_phsp.numEvt_total; i++){
        sigma+=func0(evts_phsp.costhetas[i], evts_phsp.sinthetas[i], evts_phsp.cosphis[i], evts_phsp.sinphis[i], par);
    }
    return sigma;
}


double func0 (
        vector<double> costheta, vector<double> sintheta,
        vector<double> cosphi,   vector<double> sinphi,
        double *par){

    double  dGdOm = 0;

    double  alpha = par[alpha_num];
    double    phi = par[phi_num];
    double   beta = cal_beta(alpha,phi);
    double  gamma = cal_gamma(alpha,phi);
    double alpha1 = par[alpha1_num];
    double   phi1 = par[phi1_num];
    double  beta1 = cal_beta(alpha1,phi1);
    double gamma1 = cal_gamma(alpha1,phi1);
    double alpha2 = par[alpha2_num];
    double     Px = par[Px_num];
    double     Py = par[Py_num];
    double     Pz = par[Pz_num];


    double cT0 = costheta[0];
    double cT1 = costheta[1];
    double cT2 = costheta[2];
    double sT0 = sintheta[0];
    double sT1 = sintheta[1];
    double sT2 = sintheta[2];

    double cP0 = cosphi[0];
    double cP1 = cosphi[1];
    double cP2 = cosphi[2];
    double sP0 = sinphi[0];
    double sP1 = sinphi[1];
    double sP2 = sinphi[2];

    double K1 = sT0 * (Px*cP0 + Py*sP0) + Pz*cT0;
    double K2 = cT0 * (Px*cP0 + Py*sP0) - Pz*sT0;
    double K3 = Py*cP0 - Px*sP0;
 
    double F_T1 = 0.25*( 9.*pow(cT1,2) + 1 );
    double G_T1 = 1.5*cT1;
    double f    = 0.5*( 3.*cT1 - 1. );
    double g    = 0.5*( 3.*cT1 + 1. );
   
    double AplusD  =  F_T1 * ( 1 + alpha*alpha1*cT1 + K1*(alpha + alpha1*cT1) ) 
                    - G_T1 * ( alpha*alpha1   + cT1 + K1*(alpha1 + alpha*cT1) )
                    + alpha1*f*g*sT1 * ( K2*(gamma*cP1 - beta*sP1) + K3*(gamma*sP1 + beta*cP1) );

    double AminusD =  F_T1 * ( alpha1 + alpha*cT1 + K1*(alpha*alpha1 + cT1) )
                    - G_T1 * ( alpha + alpha1*cT1 + K1*(1 + alpha*alpha1*cT1) )
                    +        f*g*sT1 * ( K2*(gamma*cP1 - beta*sP1) + K3*(gamma*sP1 + beta*cP1) );

    double douReB  = gamma1 * ( 
                                - f*g*(alpha + K1)*sT1
                                + K2*    ( F_T1*(gamma*sP1 + beta*cP1) + G_T1*(gamma*cP1 - beta*sP1) )
                                + K2*cT1*( F_T1*(gamma*cP1 - beta*sP1) - G_T1*(gamma*sP1 + beta*cP1) )
                            )
                    + beta1 * (
                                + K3*    ( F_T1*(gamma*cP1 - beta*sP1) - G_T1*(gamma*sP1 + beta*cP1) )
                                + K3*cT1*( F_T1*(gamma*sP1 + beta*cP1) + G_T1*(gamma*cP1 - beta*sP1) )
                            );
    
    double douImB  = beta1 * (
                                - f*g*(alpha + K1)*sT1
                                + K2*    ( F_T1*(gamma*sP1 + beta*cP1) + G_T1*(gamma*cP1 - beta*sP1) )
                                + K2*cT1*( F_T1*(gamma*cP1 - beta*sP1) - G_T1*(gamma*sP1 + beta*cP1) )
                            )
                    - gamma1 * (
                                + K3*    ( F_T1*(gamma*cP1 - beta*sP1) - G_T1*(gamma*sP1 + beta*cP1) )
                                + K3*cT1*( F_T1*(gamma*sP1 + beta*cP1) + G_T1*(gamma*cP1 - beta*sP1) )
                            );
    
    dGdOm = AplusD + alpha2*AminusD*cT2 + alpha2*( douReB*cP2 - douImB*sP2 )*sT2;
    

    return dGdOm;
}



Event readData(double alpha,double Phi,double alpha1,double Phi1,double alpha2,double Px,double Py,double Pz,int total_events,long SEED){

RooRandom::randomGenerator()->SetSeed(SEED);
double purity=0.75;//fraction of signal events and total events
double efficiency=1;
total_events=total_events*efficiency;
total_events=gRandom->Poisson(total_events);//total event Poison distribution

//initial value

double beta_Omegac=sqrt(1-pow(alpha,2))*sin(Phi/180*TMath::Pi());
double beta_Omegam=sqrt(1-pow(alpha1,2))*sin(Phi1/180*TMath::Pi());

double gamma_Omegac=sqrt(1-pow(alpha,2))*cos(Phi/180*TMath::Pi());
double gamma_Omegam=sqrt(1-pow(alpha1,2))*cos(Phi1/180*TMath::Pi());


RooRealVar Phi_Omegac("Phi_Omegac","Phi_Omegac",Phi,-180,180);

RooRealVar Phi_Omegam("Phi_Omegam","Phi_Omegam",Phi1,-180,180);

RooRealVar alphaOmegac("alphaOmegac","alphaOmegac",alpha,-1,1);
RooFormulaVar betaOmegac("betaOmegac","sqrt(1-@0*@0)*sin(@1*TMath::Pi()/180)",RooArgList(alphaOmegac,Phi_Omegac));
RooFormulaVar gammaOmegac("gammaOmegac","sqrt(1-@0*@0)*cos(@1*TMath::Pi()/180)",RooArgList(alphaOmegac,Phi_Omegac));

RooRealVar alphaOmegam("alphaOmegam","alphaOmegam",alpha1);//+- 0.010
RooFormulaVar betaOmegam("betaOmegam","sqrt(1-@0*@0)*sin(@1*TMath::Pi()/180)",RooArgList(alphaOmegam,Phi_Omegam));
RooFormulaVar gammaOmegam("gammaOmegam","sqrt(1-@0*@0)*cos(@1*TMath::Pi()/180)",RooArgList(alphaOmegam,Phi_Omegam));

RooRealVar alphaLam("alphaLam","alphaLam",alpha2);

RooRealVar P_x("P_x","P_x",Px,-1,1);
RooRealVar P_y("P_y","P_y",Py,-1,1);
RooRealVar P_z("P_z","P_z",Pz,-1,1);

RooRealVar theta0("theta0","theta0",0,TMath::Pi());
RooRealVar theta1("theta1","theta1",0,TMath::Pi());
RooRealVar theta2("theta2","theta2",0,TMath::Pi());
RooRealVar phi0("phi0","phi0",-1*TMath::Pi(),TMath::Pi());
RooRealVar phi1("phi1","phi1",-1*TMath::Pi(),TMath::Pi());
RooRealVar phi2("phi2","phi2",-1*TMath::Pi(),TMath::Pi());



 RooAngularPdf *sigpdf= new RooAngularPdf("sigpdf","angluar", phi0, theta0, phi1, theta1, phi2, theta2, alphaOmegac,  betaOmegac,  gammaOmegac, alphaOmegam, betaOmegam, gammaOmegam, alphaLam,  P_x, P_y, P_z);
RooDataSet *AngData = sigpdf->generate(RooArgSet( phi0, theta0, phi1, theta1, phi2, theta2), total_events);

TTree *t_data=AngData->GetClonedTree();


    int total = t_data->GetEntries();

    Double_t var1, var2, var3, var4, var5, var6;
    t_data->SetBranchAddress("theta0", &var1);
    t_data->SetBranchAddress("theta1", &var2);
    t_data->SetBranchAddress("theta2", &var3);
    t_data->SetBranchAddress("phi0",   &var4);
    t_data->SetBranchAddress("phi1",   &var5);
    t_data->SetBranchAddress("phi2",   &var6);


    vector<vector<double> >thetas, phis; thetas.clear(); phis.clear();
    vector<double> theta, phi;
    
    for(int ii=0; ii<total; ii++){
        t_data->GetEntry(ii);

        theta.clear();
        theta.push_back(var1);
        theta.push_back(var2);
        theta.push_back(var3);
        thetas.push_back(theta);

        phi.clear();
        phi.push_back(var4);
        phi.push_back(var5);
        phi.push_back(var6);
        phis.push_back(phi);

    }
    delete t_data;
    return Event(thetas, phis);
}

double cal_beta(double alpha,double phi)
{
    return sqrt(1. - pow( alpha ,2) )*TMath::Sin(phi*TMath::Pi()/180.);
}

double cal_gamma(double alpha,double phi)
{
    return sqrt(1. - pow( alpha ,2) )*TMath::Cos(phi*TMath::Pi()/180.);
}

TTree *save_result(TTree *t_data,double *par,double *err,double *betaerr,double *gammaerr,int flag)
{
            double  alpha = par[alpha_num];
            double    phi = par[phi_num];
            double alpha1 = par[alpha1_num];
            double   phi1 = par[phi1_num];
            double alpha2 = par[alpha2_num];
            double     Px = par[Px_num];
            double     Py = par[Py_num];
            double     Pz = par[Pz_num];

            double  alpha_err = err[alpha_num];
            double    phi_err = err[phi_num];
            double alpha1_err = err[alpha1_num];
            double   phi1_err = err[phi1_num];
            double alpha2_err = err[alpha2_num];
            double     Px_err = err[Px_num];
            double     Py_err = err[Py_num];
            double     Pz_err = err[Pz_num];

            double   beta = cal_beta(alpha,phi);
            double  gamma = cal_gamma(alpha,phi);
            double  beta1 = cal_beta(alpha1,phi1);
            double gamma1 = cal_gamma(alpha1,phi1);

            double  beta_err=betaerr[0];
            double gamma_err=gammaerr[0];
            double  beta1_err=betaerr[1];
            double gamma1_err=gammaerr[1];


            double diff[8];
            for (int i=0;i<8;i++)
            {
                diff[i]=par[i]-gen_settings[i];

            }
            double  alpha_diff = diff[alpha_num];
            double    phi_diff = diff[phi_num];
            double alpha1_diff = diff[alpha1_num];
            double   phi1_diff = diff[phi1_num];
            double alpha2_diff = diff[alpha2_num];
            double     Px_diff = diff[Px_num];
            double     Py_diff = diff[Py_num];
            double     Pz_diff = diff[Pz_num];

        double beta_initial=cal_beta(gen_settings[alpha_num],gen_settings[phi_num]);
        double gamma_initial=cal_gamma(gen_settings[alpha_num],gen_settings[phi_num]);
        double beta1_initial=cal_beta(gen_settings[alpha1_num],gen_settings[phi1_num]);
        double gamma1_initial=cal_gamma(gen_settings[alpha1_num],gen_settings[phi1_num]);

        double beta_diff = beta-beta_initial;
        double beta1_diff = beta1-beta1_initial;
        double gamma_diff = gamma-gamma_initial;
        double gamma1_diff = gamma1-gamma1_initial;

if (flag==0)
{
            t_data->Branch("Omegac_alpha",&alpha,"Omegac_alpha/D");
            t_data->Branch("Omegac_alpha_diff",&alpha_diff,"Omegac_alpha_diff/D");
            t_data->Branch("Omegac_alphaError",&alpha_err,"Omegac_alphaError/D");

            t_data->Branch("Omegac_Phi",&phi,"Omegac_Phi/D");
            t_data->Branch("Omegac_PhiError",&phi_err,"Omegac_PhiError/D");
            t_data->Branch("Omegac_Phi_diff",&phi_diff,"Omegac_Phi_diff/D");

            t_data->Branch("Omegac_beta",&beta,"Omegac_beta/D");
            t_data->Branch("Omegac_beta_diff",&beta_diff,"Omegac_beta_diff/D");
            t_data->Branch("Omegac_gamma",&gamma,"Omegac_gamma/D");
            t_data->Branch("Omegac_gamma_diff",&gamma_diff,"Omegac_gamma_diff/D");
            t_data->Branch("Omegac_betaError",&beta_err,"Omegac_betaError/D");
            t_data->Branch("Omegac_gammaError",&gamma_err,"Omegac_gammaError/D");


            t_data->Branch("Omegam_Phi",&phi1,"Omegam_Phi/D");
            t_data->Branch("Omegam_PhiError",&phi1_err,"Omegam_PhiError/D");
            t_data->Branch("Omegam_Phi_diff",&phi1_diff,"Omegam_Phi_diff/D");

            t_data->Branch("Omegam_beta",&beta1,"Omegam_beta/D");
            t_data->Branch("Omegam_beta_diff",&beta1_diff,"Omegam_beta_diff/D");
            t_data->Branch("Omegam_gamma",&gamma1,"Omegam_gamma/D");
            t_data->Branch("Omegam_gamma_diff",&gamma1_diff,"Omegam_gamma_diff/D");
            t_data->Branch("Omegam_betaError",&beta1_err,"Omegam_betaError/D");
            t_data->Branch("Omegam_gammaError",&gamma1_err,"Omegam_gammaError/D");


            t_data->Branch("P_x",&Px,"P_x/D");
            t_data->Branch("P_y",&Py,"P_y/D");
            t_data->Branch("P_z",&Pz,"P_z/D");

            t_data->Branch("P_xError",&Px_err,"P_xError/D");
            t_data->Branch("P_yError",&Py_err,"P_yError/D");
            t_data->Branch("P_zError",&Pz_err,"P_zError/D");

            t_data->Branch("P_x_diff",&Px_diff,"P_x_diff/D");
            t_data->Branch("P_y_diff",&Py_diff,"P_y_diff/D");
            t_data->Branch("P_z_diff",&Pz_diff,"P_z_diff/D");
}

            t_data->Fill();
return t_data;



}

